Immunoinformatics Study: Multi-Epitope Based Vaccine Design from SARS-CoV-2 Spike Glycoprotein

The coronavirus disease 2019 outbreak has become a huge challenge in the human sector for the past two years. The coronavirus is capable of mutating at a higher rate than other viruses. Thus, an approach for creating an effective vaccine is still needed to induce antibodies against multiple variants with lower side effects. Currently, there is a lack of research on designing a multiepitope of the COVID-19 spike protein for the Indonesian population with comprehensive immunoinformatic analysis. Therefore, this study aimed to design a multiepitope-based vaccine for the Indonesian population using an immunoinformatic approach. This study was conducted using the SARS-CoV-2 spike glycoprotein sequences from Indonesia that were retrieved from the GISAID database. Three SARS-CoV-2 sequences, with IDs of EIJK-61453, UGM0002, and B.1.1.7 were selected. The CD8+ cytotoxic T-cell lymphocyte (CTL) epitope, CD4+ helper T lymphocyte (HTL) epitope, B-cell epitope, and IFN-γ production were predicted. After modeling the vaccines, molecular docking, molecular dynamics, in silico immune simulations, and plasmid vector design were performed. The designed vaccine is antigenic, non-allergenic, non-toxic, capable of inducing IFN-γ with a population reach of 86.29% in Indonesia, and has good stability during molecular dynamics and immune simulation. Hence, this vaccine model is recommended to be investigated for further study.


Introduction
The COVID-19 pandemic, caused by the SARS-CoV-2 virus, has become a world-wide issue. SARS-CoV-2 is a beta coronavirus in the Coronaviridae family and is believed to have originated from bats. It is most closely related to BatCoV-RaTG13, with a sequence similarity rate of approximately 96%, and to SARS-CoV, with a similarity rate of roughly 80% [1,2]. The SARS-CoV-2 virus has a positive, and single-stranded RNA with 29.8 kilobases (kb) long. This genome encodes non-structural proteins (Open Reading Frames (ORFs) a and b) and structural proteins (Envelope (E), Spike (S), Nucleocapsid (N), and Membrane (M)) [3,4]. The spike glycoprotein of SARS-CoV-2 is comprised of two subunits: S1 (residues 14 to 685) and S2 (residues 686 to 1273) and has a total length of 1273 amino acids (AA). The S1 subunit contains N-terminal and receptor-binding domains (RBD), which play a critical role in the infection process as it binds to the ACE2 receptor in human cells [5,6]. Research has revealed the potential of the spike glycoprotein as an antigenic region [7]. This discovery makes it a promising focus for research and vaccine development. Vaccines trigger the production of specific antibodies against a disease, but the development of conventional vaccines made from the whole virus can lead to side effects [8], and potential reverse virulence [9,10]. The presence of viral mutations can also change the infection mechanism so that a new vaccine is needed to induce the specific antibodies [9,11]. In Indonesia, there are several types of COVID-19 vaccines available, including mRNA-based vaccines (Pfizer and Moderna), adenovirus vector-based (AstraZeneca), and inactivated virus vaccine (Sinovac). The main challenge in the development of the COVID-19 vaccine is the high mutation rate and transmissibility of the virus [12]. Thus, a platform that can boost vaccine design in a shorter time is crucial to minimize the adverse effect on public health conditions. To address this, immunoinformatics is being used to design vaccines more efficiently, providing a precise, rapid, and effective approach to vaccine development [13]. Utilizing large immunoinformatics databases can support multiepitope vaccine design, which has been shown to induce an immune response against a specific virus without side effects [14]. The use of peptide series in multiepitope vaccines can prevent and treat viral infections by inducing an immune response from CTL, Th, and B cells [15]. Adding a linker adjuvant between each epitope can also enhance the vaccine's immunogenic effect [10]. Therefore, an immunoinformatic study to see the potential of a multiepitope vaccine for SARS-CoV-2 can help to foresee alternatives in developing vaccines for the global population.
In this study, we performed an immunoinformatic analysis to design a multiepitope vaccine for three SARS-CoV-2 S proteins detected in Indonesia. A previous study by Gustiananda et al. (2021) [16] revealed polyprotein sequences for vaccine candidates based on the COVID-19 ORF1ab. They were validated using the physicochemical structure and analyzed further through cross-reactivity, molecular docking, and immune simulation. Meanwhile, Febrianti and Narulita (2022) [17] designed a multiepitope COVID-19 vaccine using spike glycoprotein and validated its sequence using physicochemical structure, which also continued to the design of a plasmid vector. Regarding these studies, our research started by constructing a phylogenetic tree from the evolution of the target sequence in Indonesia, searching the epitope of T and B cells, constructing and validating the vaccines and their 3D structure, docking them with TLR-3 or TLR-4, performing a molecular dynamics study, in silico immune simulations, and constructing a plasmid vector.

Study Design
The study started by retrieving Indonesia's sequence of SARS-CoV-2 spike glycoprotein from the database. These data were converted into amino acid sequences, then phylogenetic tree was constructed Three sequences were selected for further analysis against the original sequence (Wuhan-Hu-1). These sequences were used to predict CTL epitopes, HTL epitopes, B-cell epitopes, and IFN-γ. The antigenicity, allergenicity, toxicity, and population coverage were analyzed. A multiepitope vaccine was then designed, modeled, and validated. Finally, molecular docking, molecular dynamics analysis, and in silico immune simulations were performed. Figure 1 shows the design of this study.
Vaccines 2023, 11, x FOR PEER REVIEW 3 Figure 1. Epitope grouping sequence on the vaccine structure and its linkers.

The Sequences of SARS-CoV-2 Spike Glycoprotein in Indonesia
The nucleotide sequences of SARS-CoV-2 were retrieved from the Global Init on Sharing Avian Influenza Database (GISAID) (https://gisaid.org/, accessed on 13 ary 2021) [18], on 19 January 2021. The sequences obtained were the viral sequence appeared in Indonesia. Meanwhile, the nucleotide sequences of SARS-CoV-2 Wuhan 1 were retrieved from the National Centre for Biotechnology Information (N (https://www.ncbi.nlm.nih.gov/, accessed on 16 January 2021) (Accession ID: NC_04

The Sequences of SARS-CoV-2 Spike Glycoprotein in Indonesia
The nucleotide sequences of SARS-CoV-2 were retrieved from the Global Initiative on Sharing Avian Influenza Database (GISAID) (https://gisaid.org/, accessed on 13 Jan-Vaccines 2023, 11, 399 3 of 21 uary 2021) [18], on 19 January 2021. The sequences obtained were the viral sequences that appeared in Indonesia. Meanwhile, the nucleotide sequences of SARS-CoV-2 Wuhan-Hu-1 were retrieved from the National Centre for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/, accessed on 16 January 2021) (Accession ID: NC_045512). During this research, a new variant virus occurred in Indonesia in March 2021, the SARS-CoV-2 B.1.1.7. Its sequence was retrieved from GISAID and included in the analysis. All the nucleotide sequences retrieved were converted into amino acids using the MEGA-X software package [19]. Then, multiple sequence alignment was performed using ClustalW [20]. The phylogenetic tree was constructed with maximum likelihood and the neighbor-joining methods by MEGA-X [19][20][21][22]. The sequence at the outer clade of the phylogenetic tree was chosen to be further analyzed as it represents the most common sequence that spread in Indonesia (ID EPI_ISL_574613). Meanwhile, the sequences with the furthest clade were also selected as the most distinct sequence and showed the lowest similarity (ID EPI_ISL_576116).

IFN-γ Prediction
The IFN-γ stimulation test was performed using IFN-epitopes (https://webs.iiitd. edu.in/raghava/ifnepitope/, accessed on 28 February 2021) for selected HTL and B-cell epitopes. The parameters used were no-splitting at window length, the motif and SVM hybrid algorithm, also the prediction model of IFN-γ versus non-IFN-γ.

Multiepitope Vaccine Construction, Modeling, and Validation
The selected epitopes were utilized to design the multiepitope vaccine. Beta-defensin (NCBI Accession ID: P81534) was used as an adjuvant at the N-terminal [37]. The linker EAAAK was used to link the adjuvant to the CTL, and the linker AAY was used to link the CTL-CTL epitopes. The linker GPGPG was used to connect CTL and HTL epitopes, among HTL epitopes, HTL and B cell epitopes, and among B cell epitopes ( Figure 1). The 3D structure of the predicted vaccines was built using TrRosetta (https://yanglab.nankai. edu.cn/trRosetta/, accessed on 17 April 2021) [38]. Then, the structural accuracy of the vaccine model was evaluated using PROCHECK (http://www.ebi.ac.uk/thornton-srv/ software/PROCHECK/, accessed on 22 April 2021) [39][40][41]. The structural quality of the vaccine model was measured using ERRAT (https://saves.mbi.ucla.edu/, accessed on 23 April 2021) with a threshold > 50 (maximum value 100) [42].

Molecular Dynamics
Molecular dynamics simulation was performed using GROMACS and CHARMM36 force field for each vaccine model and its complex against TLR3 or TLR4. The simulation was performed using TIP3P as a water model and with 0.15 NaCl by adding the appropriate number of ions (sodium or chloride). Periodic boundary conditions (PBCs) were applied to the system in all the spatial directions. LINCS algorithms were used, and all hydrogen bonds were constrained. A 1.2 nm distance cutoff for the short-distance electrostatic and Van der Waals interactions was used. Particle Mesh Ewald algorithm (PME) was used to calculate the long-range electrostatic forces. The steepest descent algorithm was used to minimize the system's energy. The system was then allowed to reach an equilibrium state through the NVT ensemble by using the V-Rescale thermostat at 300 K, then through the NPT ensemble by using the Parrinello-Rahman barostat at 1 atm. A 50 ns simulation was performed for the individual vaccine model and 10 ns simulation was performed for the complex system.

In Silico Cloning and Immune Simulation
The designed vaccine was simulated using a C-ImmSim web server (https://kraken. iac.rm.cnr.it/C-IMMSIM/, accessed on 3 August 2021) [52,53]. The simulation parameters were kept default except for the time steps in 1, 84, and 170 with 1050 total simulation phases. Hence, there were three injections at the interval of four weeks [54]. In order to optimize the expression of the designed vaccine, codon optimization, and cloned vector were conducted. Codon optimization was done using the VectorBuilder (https://en.vectorbuilder.com/ tool/codon-optimization.html, accessed on 7 August 2021). After optimization, the vaccine sequence was inserted into the pET28a-EgC plasmid vector from the Addgene (https: //www.addgene.org/, accessed on 8 August 2021) with XbaI and XhoI restriction sites using Serial Cloner software.

Evolutionary Relationship of SARS-CoV-2 in Indonesia and Detection of Its Mutation Site
In this study, a phylogenetic tree was built to analyze the genetic changes of SARS-CoV-2 in Indonesia and understand how new virus variants may emerge. A total of 166 SARS-CoV-2 spike glycoprotein nucleotide sequences collected from GISAID and sampled from Indonesia were analyzed. The sequences shared a common ancestral lineage to the wild type of SARS-CoV-2, Wuhan-Hu-1. Based on the constructed phylogenetic tree (Table S1) the EIJK-61453 has the most similarity to the Wuhan-Hu-1, with 100% similarity, while the UGM0002 has the least similarity, with 99.84% similarity. These two sequences along with a later emerged sequence, B.1.1.7 were then used for further analysis. The mutation of SARS-CoV-2 has occurred over time and led to changes in its properties, this may affect the human antibodies' recognition. Several studies have reported that any mutation in the SARS-CoV-2 spike protein can affect neutralization.
Compared to Wuhan-Hu-1, there were mutations in UGM0002 and B.1.1.7 except EIJK-61453. The mutated amino acids are listed in Table 2. The deletion of H69 and V70 results in two-fold higher infectivity than the Wuhan-Hu-1 and decreases the sensitivity of antibodies to neutralize the virus. The N501 mutation into Y increases the ability of the virus to bind to ACE2. The D614 change into G that occurred in UGM0002 and B.1.1.7 also causes an increase in the infectivity of the SARS-CoV-2 to neutralizing antibodies [55][56][57][58]. However, Korber's study revealed that this does not imply that the virus is antibody-resistant [59]. In addition, the P681 mutation into R at the furin cleavage site has been shown to increase in the membrane fusion ability, leading to enhanced transmission of the virus in the body [60]. The impact of other mutations remains unknown.
The mutation on SARS-CoV-2 spike glycoprotein changes the segments targeted by antibodies. Several studies reported that these mutations can evade monoclonal antibodies [56,61] and influence polyclonal antibody recognition [62,63]. These mutations have been found to be present in the human population [61].

Identification and Selection of T-Cell Epitopes
A total of 300 identified CTL epitopes have been identified. The final list of the selected epitopes for each sequence is presented in Table 3. The epitopes were chosen based on low percentile and the IC 50 value. Antigenicity, allergenicity, as well as toxicity were considered. As a result, seven epitopes from each Wuhan-Hu-1 and EIJK-61453, six epitopes from UGM0002, and five epitopes from B.1.1.7 were shortlisted for further multiple-epitope-based vaccine design.  The selected HTL epitopes are listed in Table 4. There were five selected HTL epitopes for each Wuhan-Hu-1 and EIJK-61453, four epitopes for the UGM0002, and three epitopes for the B.1.1.7. The selection of epitopes was also based on the value of IFN-γ, a cytokine mainly produced by Natural Killer (NK) cells. IFN-γ plays a vital role in inhibiting viral replication. The presence of antiviral IFN-γ in APC can elevate the stimulation of the adaptive immune response to the infection, resulting in the generation of memory in the body to recognize, and neutralize the virus if a similar infection occurs in the future. The presence of IFN-γ stimulation also increases the antigen presentation process during T-cell priming, both in terms of efficiency, quantity, quality, and diversity of peptides loaded into the MHC-I receptor [64,65]. Therefore, the epitope selected for vaccine construction is preferably derived from an epitope capable of stimulating IFN-γ.

Identification and Selection of B-Cell Epitopes
B-cell epitopes consist of linear and discontinuous (conformational) epitopes with 5-20 AA lengths. Its selections were based on an antigenicity score higher than 0.5 and the IFN-γ stimulation ability. Table 5 shows selected linear B-lymphocyte (LBL) epitopes that meet the criteria. There were five selected LBL epitopes from each Wuhan-Hu-1, EIJK-61453, and B.1.1.7, also six selected LBL epitopes from UGM0002. Table S4 shows two discontinuous B-cell epitopes shortlisted from each sequence.

Indonesian Population Coverage from Selected Epitopes
Each epitope can bind to a specific HLA allele depending on the population coverage. Therefore, in this study, the coverage of the Indonesian population was analyzed to determine the specific epitope that might bind to the Indonesian HLA allele. The selected epitopes of Wuhan-Hu-1 and EIJK-61453 have a coverage score of 78.26%, while the selected epitope of UGM0002 and B.1.1.7 have 86.29% and 84.28% coverage scores, respectively. Tables S1-S3 show at least one HLA from Austronesian Indonesian ethnicity can bind to the selected CTL and HTL epitopes. The details of the coverage distribution of each epitope are attached in Figure S2.

Epitope Analysis in Human Peptide
The mimicry of the epitope needs to be considered in the construction of a multi-epitopebased COVID-19 vaccine as it can result in an auto-immune effect [66]. Clinical reports suggest that no T-cell and B-cell epitopes for the SARS-CoV-2 vaccine candidate have been found to be homologous or share similarities with proteins that cause auto-immune effects or increase the severity of COVID-19 disease [67][68][69]. Another study found that there were epitopes with similarities to human peptides on six or seven consecutive peptides based on the results of CTL analysis using sequences from the IEDB database. KIYSKHTPI and SPRRARSVA are among those epitopes [70]. However, no reports specifically mention the effect of this peptide similarity, but this can be considered for further in vitro/in vivo analysis.

Vaccine Construction and Its Validation
The multi-epitope vaccine was constructed using all the selected epitopes by optimizing the arrangement of multiple epitopes connected by linkers and an adjuvant. The possible epitope combination was then evaluated by Ramachandran plot, TrRosetta, and ERRAT score. Three vaccines that fulfill the criteria are presented in Table S5. During the optimization, we found that the prediction score from TrRosetta can be optimized if the sequence is arranged in such a manner that is familiar to the existing protein sequence template, thus making the algorithm easier to recognize the vaccine structure. Moreover, AA length also affects the quality of the vaccine structure. The increasing number of AA will potentially reduce the structure's quality in Psi-Phi stability and the TM-Score and ERRAT (Table S5). Table 6 displays the optimal result of the vaccine construction prediction. Wuhan-Hu-1 and EIJK-61453 were found to be identical, only the EIJK sequence will be subjected to further analysis. The 3D visualization and Ramachandran plot of those are shown in Figures 2 and 3, respectively. The use of linkers in multi-epitope vaccine development can primarily minimize the probability of deformation in the vaccine structure [71]. We utilized EAAAK, AAY, and GPGPG sequence linkers to build the complex vaccine and determine the best quality of the produced vaccine. The use of different linkers may also impart a specific role to the vaccine structure. For instance, the EAAAK linker avoids interference from other proteins when an adjuvant binds to its receptor. The other linkers such as AY and GPGPG were used to induce the immune response and provide a site for proteasomal cleavage. Furthermore, glycine and proline are often present in a stable beta turn. Proline's cyclic structure is ideal for beta-turn, while glycine has the smallest side chain compared to all other AAs making it the most sterically flexible linker [72].
We also use an adjuvant to enhance the efficacy of our designed vaccine. In comparison to inactivated or attenuated virus-based vaccine, a peptide-based vaccine has lower immunogenicity, and therefore an adjuvant is needed. It can boost the immune response and acts as a vehicle to deliver the vaccine structure [10]. Moreover, different adjuvants can trigger varying immune responses. The Matrix-M™ adjuvant elicits TH1 immune response, while the AS03 adjuvant elicits both Th1 and Th2 cytokine responses [10]. In our study, Beta-defensin was used as an adjuvant and functioned as an immunomodulator and anti-microbial agent [73]. possible epitope combination was then evaluated by Ramachandran plot, TrRosetta, and ERRAT score. Three vaccines that fulfill the criteria are presented in Table S5. During the optimization, we found that the prediction score from TrRosetta can be optimized if the sequence is arranged in such a manner that is familiar to the existing protein sequence template, thus making the algorithm easier to recognize the vaccine structure. Moreover, AA length also affects the quality of the vaccine structure. The increasing number of AA will potentially reduce the structure's quality in Psi-Phi stability and the TM-Score and ERRAT (Table S5). Table 6 displays the optimal result of the vaccine construction prediction. Wuhan-Hu-1 and EIJK-61453 were found to be identical, only the EIJK sequence will be subjected to further analysis. The 3D visualization and Ramachandran plot of those are shown in Figures 2 and 3, respectively. The use of linkers in multi-epitope vaccine development can primarily minimize the probability of deformation in the vaccine structure [71]. We utilized EAAAK, AAY, and GPGPG sequence linkers to build the complex vaccine and determine the best quality of the produced vaccine. The use of different linkers may also impart a specific role to the vaccine structure. For instance, the EAAAK linker avoids interference from other proteins when an adjuvant binds to its receptor. The other linkers such as AY and GPGPG were used to induce the immune response and provide a site for proteasomal cleavage. Furthermore, glycine and proline are often present in a stable beta turn. Proline's cyclic structure is ideal for beta-turn, while glycine has the smallest side chain compared to all other AAs making it the most sterically flexible linker [72].
We also use an adjuvant to enhance the efficacy of our designed vaccine. In comparison to inactivated or attenuated virus-based vaccine, a peptide-based vaccine has lower immunogenicity, and therefore an adjuvant is needed. It can boost the immune response and acts as a vehicle to deliver the vaccine structure [10]. Moreover, different adjuvants can trigger varying immune responses. The Matrix-M™ adjuvant elicits TH1 immune response, while the AS03 adjuvant elicits both Th1 and Th2 cytokine responses [10]. In our study, Beta-defensin was used as an adjuvant and functioned as an immunomodulator and anti-microbial agent [73]. Vaccines 2023, 11, x FOR PEER REVIEW 11 of 23 possible epitope combination was then evaluated by Ramachandran plot, TrRosetta, and ERRAT score. Three vaccines that fulfill the criteria are presented in Table S5. During the optimization, we found that the prediction score from TrRosetta can be optimized if the sequence is arranged in such a manner that is familiar to the existing protein sequence template, thus making the algorithm easier to recognize the vaccine structure. Moreover, AA length also affects the quality of the vaccine structure. The increasing number of AA will potentially reduce the structure's quality in Psi-Phi stability and the TM-Score and ERRAT (Table S5). Table 6 displays the optimal result of the vaccine construction prediction. Wuhan-Hu-1 and EIJK-61453 were found to be identical, only the EIJK sequence will be subjected to further analysis. The 3D visualization and Ramachandran plot of those are shown in Figures 2 and 3, respectively. The use of linkers in multi-epitope vaccine development can primarily minimize the probability of deformation in the vaccine structure [71]. We utilized EAAAK, AAY, and GPGPG sequence linkers to build the complex vaccine and determine the best quality of the produced vaccine. The use of different linkers may also impart a specific role to the vaccine structure. For instance, the EAAAK linker avoids interference from other proteins when an adjuvant binds to its receptor. The other linkers such as AY and GPGPG were used to induce the immune response and provide a site for proteasomal cleavage. Furthermore, glycine and proline are often present in a stable beta turn. Proline's cyclic structure is ideal for beta-turn, while glycine has the smallest side chain compared to all other AAs making it the most sterically flexible linker [72].
We also use an adjuvant to enhance the efficacy of our designed vaccine. In comparison to inactivated or attenuated virus-based vaccine, a peptide-based vaccine has lower immunogenicity, and therefore an adjuvant is needed. It can boost the immune response and acts as a vehicle to deliver the vaccine structure [10]. Moreover, different adjuvants can trigger varying immune responses. The Matrix-M™ adjuvant elicits TH1 immune response, while the AS03 adjuvant elicits both Th1 and Th2 cytokine responses [10]. In our study, Beta-defensin was used as an adjuvant and functioned as an immunomodulator and anti-microbial agent [73].

Docking Study of the Constructed Vaccine Model against TLR-3
All three vaccine models were successfully docked with TLR-3, as shown in Figure 4. The details of the binding interaction are shown in Table 7. The binding affinity indicates the binding interaction of the two molecules. The docked complex of Wuhan-Hu-1-TLR-3 and EIJK-61453-TLR-3 present the exact same model with the binding affinity of −18.7 kcal/mol, and 0.5 ± 0.3 RMSD score. Meanwhile, the complex of UGM0002-TLR-3 and B.

Docking Study of the Constructed Vaccine Model against TLR-3
All three vaccine models were successfully docked with TLR-3, as shown in Figure  4. The details of the binding interaction are shown in Table 7  TLR is a type of Pattern Recognition Receptor (PRRs) that functions as an early detector of pathogens. Signaling in the TLR will induce an innate and adaptive immune  Table 7. The bond type and amino acid residues interacting between vaccine models and TLR-3.

Vaccine Model Binding Affinity ∆G (kcal/mol) RMSD Interaction with TLR-3
EIJK-61453 Total Interactions: TLR is a type of Pattern Recognition Receptor (PRRs) that functions as an early detector of pathogens. Signaling in the TLR will induce an innate and adaptive immune response against specific antigens. TLR-3 can recognize the presence of viral infection, therefore the existence of a strong complex between the vaccine model and TLR-3 is expected. Amino acid residues GLY360, GLY361, SER362, THR363, LEU364, GLU368, HIS539, ASN541, GLU557, THR559, SER562, GLU564, GLU565, and SER566 in TLR-3 can stimulate TLR-3 signaling. Among the residues, ASN541 facilitates the TLR-3 recognition to its target, while HIS539 plays a role in TLR-3 activation.
Among all the docked models, only the EIJK-61453 vaccine-TLR-3 complex has interactions with ASN541 (neutral), through PHE12 (aromatic) with non-bonded contacts. This interaction may occur between the CG atoms, OD1, and ND2, with an average distance of 3.57Å. Hydrogen bonds also occurred between TLR-3 HIS539 (positive) and LEU15 (aliphatic) with a distance of 2.70Å. The EIJK-61453 vaccine model also interacts with TLR-3 ILE566 at residues GLY137 and GLY138 in the form of nonbonded contacts. The interaction with residue ILE566 triggers the interaction of TICAM-1 and the transduction of IFN-β and NF-kB as an innate immune response. However, the docked complex UGM0002-TLR-3 and B.1.1.7-TLR-3 have no residues bound to TLR-3 HIS539 or ASN541.
The EIJK-61453-TLR-3 complex model has the most hydrogen bonds (18 hydrogen bonds) compared to the B.1.1.7-TLR-3 (11 hydrogen bonds) and UGM0002-TLR-3 model (12 hydrogen bonds). The B.1.1.7-TLR-3 model has four salt bridge interactions while the EIJK-61453-TLR-3 and UGM0002-TLR-3 have three salt bridge interactions. Intra-protein interactions in these complexes are different from the surface interactions between proteins. In intra-protein interactions, the increased number of hydrogen bonds contribute to the stability of the protein structure, while salt bridges can enhance the folding stability of the β-helix protein structure [74][75][76][77]. This bond's degree of conformational freedom in interactions between protein surfaces is more rigid, although not as large as intra-protein interactions. Therefore, the analysis of these bindings is more inclined towards specificity for the binding site than TLR-3. Nevertheless, hydrogen bonds and salt bridges contribute to the stability of protein-protein interactions, but they do not bind specifically to the TLR-3 binding site. This interaction helps strengthen the bond between the spike glycoprotein surfaces, making it more difficult to separate [78].

Docking Study of the Constructed Vaccine Model against TLR-4
All three vaccine models were successfully docked with TLR-4. The docked complex of the EIJK-61453-TLR-4 model had the most hydrogen bonds (15 hydrogen bonds) compared to the UGM0002-TLR-4 (13 hydrogen bonds) and B.1.1.7-TLR-4 models (12 hydrogen bonds). The complex of EIJK-61453-TLR-4, UGM0002-TLR-4, and B.1.1.7-TLR-4 also had two, three, and four salt bridge bonds, respectively. There was no significant difference among the docked complex vaccine-TLR-4 in terms of binding affinity, which was in the range of −16.3 to −15.5. Figure 5 shows the docking visualization of (a) EIJK-61453 (b) UGM0002 (c) B1.1.7 vaccine model against TLR-4. The overall interactions among vaccines model-TLR-4 are presented in Table 8. The EIJK-61453-TLR-3 complex model has the most hydrogen bonds (18 hydrogen bonds) compared to the B.1.1.7-TLR-3 (11 hydrogen bonds) and UGM0002-TLR-3 mode (12 hydrogen bonds). The B.1.1.7-TLR-3 model has four salt bridge interactions while th EIJK-61453-TLR-3 and UGM0002-TLR-3 have three salt bridge interactions. Intra-protein interactions in these complexes are different from the surface interactions between pro teins. In intra-protein interactions, the increased number of hydrogen bonds contribute t the stability of the protein structure, while salt bridges can enhance the folding stability of the β-helix protein structure [74][75][76][77]. This bond's degree of conformational freedom in interactions between protein surfaces is more rigid, although not as large as intra-protein interactions. Therefore, the analysis of these bindings is more inclined towards specificity for the binding site than TLR-3. Nevertheless, hydrogen bonds and salt bridges contribut to the stability of protein-protein interactions, but they do not bind specifically to the TLR 3 binding site. This interaction helps strengthen the bond between the spike glycoprotein surfaces, making it more difficult to separate [78].

Docking Study of the Constructed Vaccine Model against TLR-4
All three vaccine models were successfully docked with TLR-4. The docked comple of the EIJK-61453-TLR-4 model had the most hydrogen bonds (15 hydrogen bonds) com pared to the UGM0002-TLR-4 (13 hydrogen bonds) and B.1.1.7-TLR-4 models (12 hydro gen bonds). The complex of EIJK-61453-TLR-4, UGM0002-TLR-4, and B.1.1.7-TLR-4 also had two, three, and four salt bridge bonds, respectively. There was no significant differ ence among the docked complex vaccine ̶ TLR-4 in terms of binding affinity, which was in the range of −16.3 to −15.5. Figure 5 shows the docking visualization of (a) EIJK-61453 (b UGM0002 (c) B1.1.7 vaccine model against TLR-4. The overall interactions among vaccine model ̶ TLR-4 are presented in Table 8.  Table 8. The bond type and amino acid residues interacting between vaccine models and TLR-4.

Vaccine Model
Binding Affinity ΔG RMSD Interaction with TLR-4

Molecular Dynamics
Molecular dynamics analysis was performed to evaluate the stability of each vaccine-TLR complex model. Figure 6 shows the RMSD of the complexes in the 10 ns molecular dynamics simulation. As shown in Figure 6a, the UGM002 and EIJK vaccines showed an earlier stability time than B117 when interacting with TLR-3, although there were some fluctuations in the EIJK vaccine from 8 to 10 ns. In contrast, Figure 6b showed that the EIJK vaccine became stable at a later time than UGM002 and B117 vaccines when interacting with TLR-4. Furthermore, Figure 7 shows the RMSF of vaccine-TLR complex model. The interaction patterns with TLR-3 and TLR-4 were almost similar for all vaccine models. There were some fluctuations from the initial residue which lasted around the 100th residue. However, among the vaccine models, the B117 vaccine display the highest level of fluctuation.

Molecular Dynamics
Molecular dynamics analysis was performed to evaluate the stability of each vaccine-TLR complex model. Figure 6 shows the RMSD of the complexes in the 10 ns molecular dynamics simulation. As shown in Figure 6a, the UGM002 and EIJK vaccines showed an earlier stability time than B117 when interacting with TLR-3, although there were some fluctuations in the EIJK vaccine from 8 to 10 ns. In contrast, Figure 6b showed that the EIJK vaccine became stable at a later time than UGM002 and B117 vaccines when interacting with TLR-4. Furthermore, Figure 7 shows the RMSF of vaccine-TLR complex model. The interaction patterns with TLR-3 and TLR-4 were almost similar for all vaccine models. There were some fluctuations from the initial residue which lasted around the 100th residue. However, among the vaccine models, the B117 vaccine display the highest level of fluctuation.

In Silico Cloning and Immune Simulations
From the immune simulation result (Figures 8-10), three candidate vaccines (EIJK-61453, UGM0002, and B.1.1.7) showed similar patterns regarding the induction of immune cells (Figures 8-10). The EIJK-61453 vaccine showed slightly higher IgM and IgG antibody responses after the second and tertiary injections than the other candidates, while the UGM0002 has the highest active B cell production. The immune responses from candidate vaccines were generally constructed from the initial injection and resulted in higher immune responses in the second and tertiary injections. This was followed by higher B cell and Th cell populations. The main distinct difference between the initial and the subsequent injections was the total population of memory cells. From the picture, it was clear that memory cells existence is higher after both secondary and tertiary injections than in the first one. Al Zamane et al. (2021) [79] stated that the next injection after the initial phase will activate B cells and T cells and they will create long-lasting protection and memory formation followed by antigen clearance during exposures. Furthermore, from the immune simulation result, innate and cytokine based immune regulation were also involved. Macrophages as an innate immune system are responsible for creating an inflammatory response to antigen exposure. Thus, its activation will release some cytokines such as IL-12 and IL-18 to stimulate adaptive immune response and IFN-gamma production to enhance immune activation [80,81]. Therefore, these candidates are capable of inducing good immune responses.

In Silico Cloning and Immune Simulations
From the immune simulation result (Figures 8-10), three candidate vaccines (EIJK-61453, UGM0002, and B.1.1.7) showed similar patterns regarding the induction of immune cells (Figures 8-10). The EIJK-61453 vaccine showed slightly higher IgM and IgG antibody responses after the second and tertiary injections than the other candidates, while the UGM0002 has the highest active B cell production. The immune responses from candidate vaccines were generally constructed from the initial injection and resulted in higher immune responses in the second and tertiary injections. This was followed by higher B cell and Th cell populations. The main distinct difference between the initial and the subsequent injections was the total population of memory cells. From the picture, it was clear that memory cells existence is higher after both secondary and tertiary injections than in the first one. Al Zamane et al. (2021) [79] stated that the next injection after the initial phase will activate B cells and T cells and they will create long-lasting protection and memory formation followed by antigen clearance during exposures. Furthermore, from the immune simulation result, innate and cytokine based immune regulation were also involved. Macrophages as an innate immune system are responsible for creating an inflammatory response to antigen exposure. Thus, its activation will release some cytokines such as IL-12 and IL-18 to stimulate adaptive immune response and IFN-gamma production to enhance immune activation [80,81]. Therefore, these candidates are capable of inducing good immune responses.     Alongside immune simulation, codon optimization was performed for the three vaccine candidates. In this work, the codon optimization used was based on Escherichia coli strain K12 expression systems. E. coli expression systems are well-known systems and still dominate in recombinant protein production. The use of E. coli as the recombinant vector has several advantages, including fast growth, ease of achieving high cell density, and the availability of components for rich complex media [82,83]. The important parameters for codon optimization are GC contents (30-70%), and a codon adaptation index (CAI) (>0.8). Table S6 showed that the %GC contents for all candidates were around 59% with CAI Alongside immune simulation, codon optimization was performed for the three vaccine candidates. In this work, the codon optimization used was based on Escherichia coli strain K12 expression systems. E. coli expression systems are well-known systems and still dominate in recombinant protein production. The use of E. coli as the recombinant vector has several advantages, including fast growth, ease of achieving high cell density, and the availability of components for rich complex media [82,83]. The important parameters for codon optimization are GC contents (30-70%), and a codon adaptation index (CAI) (>0.8). Table S6 showed that the %GC contents for all candidates were around 59% with CAI above 0.95. CAI has been used for the assessment of highly expressed genes, the closer the value of CAI to 1, the higher the chance of the genes being expressed [84]. From the results, %GC contents and CAI from vaccine candidates represent the ideal sequence for gene expression. To simulate vaccine expression, plasmid cloning vector pET28a-EgC was constructed with XbaI and XhoI as restriction sites. The plasmid constructions result revealed that EIJK-61453, UGM0002, and B.1.1.7 have 5636, 5652, and 5582 base pairs, respectively ( Figure 11).

Conclusions
Our study has revealed a stable vaccine design that was created using the SARS-CoV-2 spike glycoprotein from Indonesia, which also became our limitation of the study. The evolution of the SARS-CoV-2 virus in Indonesia is not significantly different from its wild type, Wuhan-Hu-1, with the furthest similarity value of 99.84%. Several amino acid mu-

Conclusions
Our study has revealed a stable vaccine design that was created using the SARS-CoV-2 spike glycoprotein from Indonesia, which also became our limitation of the study. The evolution of the SARS-CoV-2 virus in Indonesia is not significantly different from its wild type, Wuhan-Hu-1, with the furthest similarity value of 99.84%. Several amino acid mutations occur and affect the construction of vaccine models. In this study, the sequences of EIJK-61453, UGM0002, and B.1.1.7 were chosen and used to search for CTL, HTL, and LBL epitopes. The selected epitopes from each sequence were then arranged with linkers and an adjuvant to construct a vaccine model. Among all the vaccine models, the UGM0002 vaccine is antigenic, non-allergenic, non-toxic, capable of inducing IFN-γ with a population reach of 86.29% in Indonesia, and activates the highest rate of B cells according to immune simulations. Furthermore, molecular dynamics analysis revealed that the interactions of the vaccines with TLR-3 or TLR-4 were stable. Therefore, the vaccine model of UGM0002 is recommended for further investigations.  Table S1: Wuhan-Hu-1 and EIJK-61453 epitopes population coverage;